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Abstract 

Background: In teleosts such as Atlantic salmon {Salmo salar L), segmentation and subsequent mineralisation of 
the notochord during embryonic stages are essential for normal vertebrae formation. However, the molecular 
mechanisms leading to segmentation and mineralisation of the notochord are poorly understood. The aim of this 
study was to identif/ genes/pathways acting in gradients over time and along the anterior-posterior axis during 
notochord segmentation and immediately prior to mineralisation of the vertebral bodies in Atlantic salmon. 

Results: Notochord samples were collected from unsegmented, pre-segmented and segmented developmental 
stages. In each stage, the cellular core of the notochord was cut into three pieces along the longitudinal axis 
(anterior, mid, posterior). RNA was sequenced (22 million pair-end 100 bp/ library) and mapped to the salmon 
genome. 66569 transcripts were predicted and 55775 were annotated. In order to identif/ possible gradients leading 
to segmentation of the notochord, all 71 notochord-expressed hox genes were investigated, most of them displaying 
a typical anterior-posterior expression pattern along the notochord axis. The clustering of hox genes revealed a pattern 
that could be related to notochord segmentation. We further investigated how mineralisation is initiated in the 
notochord, and several factors related to chondrogenic lineage were identified (sox9, sox5, sox6, tgfbS, ihhb and col2al), 
suggesting a cartilage-like character of the notochord. KEGG analysis of differentially expressed genes between stages 
revealed down-regulation of pathways associated with EGM, cell division, metabolism and development at onset of 
notochord segmentation. This implies that inhibitory signals produce segmentation of the notochord. One such 
potential inhibitory signal was identified, coll la2, which was detected in segments of non-mineralising notochord. 

Conclusions: An incomplete salmon genome was successfully used to analyse RNA-seq data from the cellular core 
of the Atlantic salmon notochord. In transcriptome we found; hox gene patterns possibly linked to segmentation; 
down-regulation of pathways in the notochord at onset of segmentation; segmented expression of coll la2 in 
non-mineralised segments of the notochord; and a chondroblast-like footprint in the notochord. 
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Background 

The notochord is a midline structure that appears early in 
the embryo of all vertebrates, and has several important 
functions. It provides internal hydroskeletal support, until 
this role is taken over by the vertebral column. It also 
produces secreted factors during development that 
provide position and fate signals to adjacent ectoderm, 
paraxial mesoderm and endoderm along the dorso-ventral 
axis [1-5]. In recent years, additional functions of the 
notochord have been unravelled. Segmentation and 
subsequent mineralisation of the notochord comprise the 
initial morphogenic process in formation of the vertebral 
column in Atlantic salmon (Salmo salar L.) [6] and in 
zebrafish {Danio rerio) [2,7]. Segmented mineralisation of 
the notochord determines the localization of vertebral 
bodies, and notochord mineralisation is thereby crucial 
for normal formation of the vertebrae, as the functional 
study by Willems et al. [8] has demonstrated. Furthermore, 
pathological processes that disrupt notochord mineralisa- 
tion may lead to malformation of vertebrae [9]. However, 
mechanisms initiating segmentation and mineralisation 
within the notochord are currently unknown, and further 
studies of this topic are needed. 

The salmon embryo provides a good model for studies 
of development and differentiation, due to its slow early 
development and large size, allowing pure samples of 
the cellular core of the notochord to be isolated. This 
has facilitated some detailed molecular studies of salmon 
notochord, and the discovery of notochord-specific 
molecules such as vimentin and elastin [10]. Moreover, 
the development of the notochord and the vertebral 
column of Atlantic salmon embryos and larvae have 
previously been the subjects of detailed morphological 
studies [6,9,11-13]. This has made it possible to characterize 
in detail both notochord segmentation and the subsequent 
mineralisation that nucleates vertebral development. 

As development of the salmon notochord proceeds, 
segmentation is initially observed in the notochord 
epithelium as metameric bands of cells with alternating 
cell-axis orientations [6] concurrent with expression of 
alkaline phosphatase (ALP) in every second band [12], 
both processes defining the future positions of vertebral 
and intervertebral segments. In the following stage, 
possibly in order to permit further growth of the vertebrae, 
sclerotome-derived osteoblasts differentiate on the surface 
of the notochord centrae to form bone [9,11]. 

Despite the importance of the notochord in the initial 
shaping of the vertebral column, the genes involved and the 
network of transcriptional regulators controlling them are 
still incompletely characterized. An overall transcriptome 
assessment using RNA-seq on critical stages during the 
early segmentation stage of salmon notochord might reveal 
pathways associated with the segmentation and subsequent 
mineralisation processes. Possible factors associated with 



segmentation of the notochord could be members of the 
hox gene family, as these molecules are known to confer 
patterning formation along the anterior-posterior body axis 
in all vertebrates. In zebrafish, four hox genes have been 
found to be differentially expressed along the anterior- 
posterior axis of the notochord [14], which indicates that 
these genes may play a role in metameric morphogenic 
processes within the notochord of teleosts. In Atlantic 
salmon, 118 hox genes have been identified, including 
8 pseudogenes [15], but none of them have been studied 
with regard to notochord development. Furthermore, the 
differentiation of chordocytes, and specific features of 
notochord mineralisation, might posibly be elucidated via 
a transcriptome study, as factors typical of chondrogenesis 
and osteogenesis could be identified. For instance, 
both the hedgehog and Wnt signalling pathways play 
important roles in formation of the notochord sheath 
and the patterning of the vertebral column during early 
development in the mouse [16]. 

The primary objective of the present study was to relate 
the transcriptome of the notochord to morphological events 
in un-segmented notochord (510 day), pre-segmented 
notochord (610 day) and segmented notochord to the 
ongoing mineralisation of the sheath (710 day). This study 
employed high-throughput RNA sequencing technology 
on dissected anterior-posterior notochord segments to 
explore, at global level, the pathways and genes expressed 
in the notochord that might contribute to notochord 
development, segmentation and subsequent mineralisation 
in salmon and related species. 

Methods 

Rearing of embryos and notochord dissection 

Salmon embryos were collected from a local hatchery 
(Marine Harvest, Alvoen, Norway). Yolk-sac embryos 
were kept in a refrigerator at around 4°C during the 
sampling period. Developmental stages were classified 
by day degrees (day), which are defined as the sum of 
the daily mean ambient water temperatures (°C) for each 
day of development. Based on our earlier morphologic 
observations [6], three development stages were selected 
spanning from the un-segmented notochord, through 
the segmental differentiation of notochord epithelium, 
to the pre-mineralisation stage of the notochordal sheath 
(Figure lA): 510 ± 20 d° (Tl:un-segmented notochord); 
610 ±20 d° (T2:pre-segmented notochord); 710 ± 20 d° 
(T3: segmented notochord/mineralisation in the sheath). 
Pure fractions of notochord cells were detached from 
the fibrous sheath, following the procedure described 
by Sagstad et al. [10]. At each stage, the notochord 
was dissected into three pieces along the longitudinal 
axis (A: anterior, M: mid, P: posterior) (Figure lA), 
before being fresh-frozen in liquid nitrogen and 
stored at -80°C. The notochord samples were intact 
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Figure 1 Salmon notochord sampling scheme and morphology. (A) Notochord sampling at three developmental stages calculated in day 
degrees (d°) (Tl:510±20 d°, T2: 610 ±20 d°, T3: 710 ±20 d°), and further division into segments (A: anterior, M: mid, P: posterior) along the 
longitudinal axis. (B) A freshly dissected notochord sample, from stage T3. Stereoscope photo, scale bar: 500 um. (C D): Scanning electron 
micrographs of dissected notochord. Notice the smooth, sheath-free surface and the large central chordocytes and the epithelium germ layer 

(chordoblasts) in transverse section (D). Scale bar: 400 um in (C) and 100 um in (D). 

L J 



with traces neither of the sheath nor of surrounding 
tissues (Figure 1B,C). 

Histology 

Samples of notochord with surrounding tissues for light 
microscopy, and pure notochord fractions for scanning 
electron microscopy (SEM), were fixed and processed as 
described by Grotmol et al [6]. Histological slides were 
photographed in an Olympus Vanox X3 microscope 
(Olympus, Tokyo, Japan), while samples for SEM were 
studied in a Zeiss FESEM Supra 50 field-emission 
scanning EM (Oberkochen, Germany). 

Total RNA extraction and transcriptome sequencing 

Total RNA was extracted from notochord pieces using 
Rneasy Mini Kit (Qiagen, Hilden, Germany), according 
to the manufacturer s instruction. Only RNA samples with 
A260/280 within 1.9-2.1 and A260/230 > 1.8 measured 
by Nanodrop (ND-1000, Thermo Scientific, Wilmington, 
USA), were used for pooling. In total, notochord pieces 
(anterior, mid-, posterior parts in separate tubes) from 36 
fish from Tl, 18 fish from T2 and 27 fish from T3 were 
pooled and prepared for nine total RNA libraries (at least 
6 [ig RNA/library, labeled as: TIA, TIM, TIP, T2A, T2M, 
T2P, T3A, T3M, T3P). RNA integrity was checked for each 
library (RIN value > 9, except for T3M) by Agilent 2100 
Bioanalyzer (RIN values for individual library were 
listed in Additional file 1: Table SI). RNA libraries were 
prepared by Fasteris Life Science (Geneva, Switzerland). 



Briefly cDNA libraries were prepared using random 
hexamer primers and samples were sequenced using 
2xl00bp paired-end high-throughput mRNA sequencing 
(RNA-seq) in single lane on lUumina Hiseq 2000 
(Illumina, San Diego, CA). Due to the time-consuming 
dissecting procedures involved and limited tissue mass, 
single libraries without replicates were sequenced. 

Bioinformatical analysis of RNA-seq 

The workflow of the bioinformatics analysis is summarised 
in Figure 2. Further, The raw data has been deposited and 
can be found at The Sequence Read Archive (SRA) at 
NCBI (Accession ID: SRA129427). Library quality control 
was performed by a titration run with 50 bp reads on 
HiSeq, followed by a manual blast of 50 fasta sequences 
against salmon EST databases and quality check of 100,000 
fastq reads for each nucleotide position from each library 
(Additional file 2: Figure SI A). Fewer than 10% fasta 
sequences were annotated as ribosomal genes, and more 
than 60% aligned to salmonids ESTs (data not shown). The 
distributions of nucleotides in all reads were checked, 
although the flrst 13 bp in all reads are biased, probably 
due to random hexamer priming in the Illumina cDNA 
library preparation protocol (Additional flle 2: Figure SIB), 
as previously reported [17]. To obtain gene expression 
values of each gene, sequences were mapped against a 
draft salmon genome contig assembly (Acc. No. 
AGKDOOOOOOOO.l). To acquire gene prediction for the 
draft salmon genome, Augustus gene prediction software 
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Figure 2 Bioinformatics analysis pipeline for RNA-seq. Sliort reads were cliecl<ed for quaiity score and nucieotide distribution, and reads were not 
trimmed. Reads were mapped to tine genome witli BWA. Predictions of coding sequences (exon) of tine salmon contig assembly were performed 
using Augustus algorithm. Gene annotations were done by homologous sequence inference against Swissprot, Uniref90 and KEGG databases. Total 
mapped reads to predicted exons were summarized and normalized using the longest transcript length among all libraries and total number of reads 
per library. 66569 transcripts/genes were identified in the raw gene list with more than three reads per gene. Differentially expressed gene (DEG) 
analysis was performed using NOISeq. DEG's were further clustered by k-mean clustering, KEGG pathway and Gene Ontology ID. 



[18] was trained using PAS A gene candidates [19]. The 
PASA gene candidates were obtained by mapping salmon 
ESTs to the draft salmon genome contig assembly. The 
Augustus gene software only predicted coding sequences 
without UTRs. Predicted genes were annotated with 
various resources, including Swissprot, Uniref90 and 
KEGG. RNA-seq reads of 100 bp were aligned to the 
genome assembly with the Burrows-Wheeler Alignment 
tool (BWA) [20] in paired-end mode. Since the genome 
reference was fragmented, the paired-end information 
was only checked in cases where both reads in a pair 
mapped to the same contig. Aligned reads were analyzed 
with biopython, biojava and picard APIs (application 
programming interface). 

The commented java source code based on the Picard 
API for extracting the read count per Augustus pre- 
dicted gene is available at https://marineseq.imr.no/ 
salmon/annot2013/. "In silico" gene expression was 
calculated as: (number of mapped reads*read length)/ 
(predicted coding sequence length). In cases where the 
mapped reads overlapped coding sequence and the UTR, 
only the part of the read that mapped to the coding 
sequence was counted. If a read mapped to a coding 
sequence with a 90% match, only 0.9 was added. Reads 
with an identity match to the genome below 90% were 
removed from the analysis. Normalization of expression 



values among libraries was completed by counting 
the total number of reads from each library and then 
recalculating, so that each library was given the same 
number of counts. 

Differential gene expression 

Differentially expressed genes (DEG) were selected by 
using NOISeq [21] with a threshold 0.8. Since we did 
not have true biological replicates, we either pooled 
three segments (A, M, P) of each developmental stages 
(labelled as sumTl, sumT2, sumT3), or three stages (Tl, 
T2, T3) for each segment (sumA, sumM, sumP). We ran 
NOISeq for both spatial and temporal comparison on 
pooled libraries, using library- normalized raw read count 
for each gene. K-mean clustering was used to group genes 
with similar patterns of expression. Gene Ontology [22] 
was used to pick up genes that were annotated with ECM. 
The output ECM gene list could be cross-checked for 
DEGs from NOISeq. KEGG pathways analysis [23] was 
performed by mapping the KEGG annotated DEGs from 
NOISeq to KEGG pathways as described in the KEGG 
Mapper tool. Both raw expressions of genes and DEGs as 
fold change were plotted in pathways, and the ratio of 
number of up-regulated genes/number of down-regulated 
genes, or vice versa, were calculated as a means of ranking 
up-/down-regulated pathways. 
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Quantitative real-time PGR 

Quantitative real-time PGR (qPCR) was performed using 
SYBR Green PGR Master Mix (Applied Biosystems Inc., 
Foster Gity, GA, USA). Total RNA of all libraries was 
extracted from dissected notochord, following the same 
procedure as described for transcriptome sequencing 
from another batch of salmon embryos [10]. Only 
RNA samples with A260/280: 1.8-2.0 were pooled 
from two additional batches of salmon embryos, and 
they were used for cDNA synthesis with Superscript 
VILO cDNA Synthesis Kit (Invitrogen), according to 
the manufacturers guidelines. After TURBO DNA-free kit 
(Ambion) treatment and sodium acetate precipitation, 
qPGR reactions were run in a 7900HT Fast Real-Time 
PGR System (Applied Biosystems Inc.). The conditions for 
all reactions were 50°G for 2 min, followed by 95°G for 
10 min, 40 cycles of 95°G for 15 sec, followed by 60°G for 
1 min. Genomic fragments from the draft salmon genome 
assembly (Acc. No. AGKDOOOOOOOO.l) were used to 
design primers and probes, as shown in Additional file 3: 
Table S2. Primers were designed using Primer Express v2.0 
software (Applied Biosystems Inc.) and Primer-BLAST 
[24]. Elongation factor 1 a (efla) was used as reference 
gene for its stable expression in salmon tissues [25], 
including notochord [10,13,26]. efla was expressed at 
the same level in all libraries (Gene expression data found 
at https://marineseq.imr.no/salmon/annot2013/). Three 
biological replicate cDNAs were made for gene expression 
analysis. Standard curve analysis was performed to confirm 
similar amplification efficiency in target genes and the 
reference gene with a validation step. A no-template-control 
(NTG) and an RNA sample without reverse transcription 
(-RT), were used to control for contamination of external 
and/or genomic DNA in reactions. Melting curve analyses 
were performed for each primer pair in order to confirm 
unique amplicon reaction. One-way ANOVA (p < 0.05) and 
Tukey s multiple comparison, were performed in GraphPad 
Prism v5 (GraphPad Software, San Diego, GA. USA) to 
detect DEG in summed libraries (sumA-P and sumTl-T3) 
following analysis of the qPGR data. 

In situ hybridization 

In situ hybridization was performed according to 
Krossoy et al. [27]. Primer sequences for collla2 are 
listed in Additional file 3: Table S2. 

Results 

Transcriptome features, annotation and quality 

Around 22 million 100 bp paired-end reads were 
collected from each library, out of which 70% of the 
reads were mapped to the draft salmon genome 
(Additional file 1: Table SI). Similar mapping percentages 
were achieved using another mapping tool, Bowtie2 
[28] (data not shown). We used at least three reads 



mapped/ Augustus-predicted genes to output 66569 
genes as the raw gene list, of which, 55775 genes 
were annotated using the Uniprot/Swissprot database. 
The whole gene set was additionally annotated with 
Uniref90 and KEGG. The Augustus gene prediction, 
predicted gene sequences and annotation files can be 
downloaded from https://marineseq.imr.no/salmon/ 
annot2013/. Out of all the mapped reads, 34.8% mapped 
to predicted coding exons in the draft salmon genome; 
32.1% mapped to regions up to 1000 bp downstream of 
predicted genes, and thus estimated the expression of 
untranslated regions (3' UTRs); 29% and 4.1% mapped to 
non-genic regions and intron regions, respectively 
(Figure 3A). The average predicted coding sequence 
length of predicted gene fragments was around 872 bp 
(Figure 3B). As a measure of the quality of all libraries, 
expression distribution was measured in each library, 
which were very similar throughout, with an average 
expression level of about 1 (Figure 3G). Pair-wise 
comparisons of the number of DEGs between summed 
notochord libraries (either spatially or temporally) using 
NOISeq are shown in Figure 3D. From sumA-P compari- 
son, 2470 DEGs were identified (1568 annotated and 902 
unknown genes), of which 1263 genes were up-regulated 
and 1207 were down-regulated during the sampling period. 
From sumTl-T3 comparison, 3147 DEGs were identified 
(2655 annotated and 483 unknown genes), of which 1531 
genes were up-regulated and 1616 down-regulated. 

qPCR validation of RNA-seq gene expression values 

A total of 23 genes were selected for validation of RNA-seq 
gene expression by independent quantitative RT-PGR. Of 
these, all except runx2 were expressed (Gt<30 cycles or 
deltaGt <8 between the target gene and the reference gene 
Elongation factor la), A high correlation was detected 
between RNA-seq and qPGR gene expression level as 
relative fold change to TIA notochord library. In general, 
most of the genes showed the same trend in both methods 
(Figure 4). The expression of four genes (fetub, tgfbl, 
bmp4, bmp6), involved in the TGFp/BMP signaling 
pathways (Figure 4A); four genes {lrp6y dkk2, tVlla, fos) 
involved in the Wnt signalling pathway (Figure 4B); three 
hox genes {hoxbSab, hoxclOaa, hoxdSaa) (Figure 4G); 
three Hedgehog genes {ihhby shh, twhh) (Figure 4D); 
three major transcriptional factors {sox9y runx3, atf4) in 
osteogenic and chondrogenic lineages of mesenchyme 
stem cells (Figure 4E); three genes {col2al, colllal, 
collla2) coding fibrillar collagen chains (Figure 4F) and 
two calcium-binding protein-coding genes {si 00a, sprc) 
(Figure 4G) were shown. Further statistical analysis 
(NOISeq and one-way ANOVA among summed libraries 
(sumA-P or sumTl-T3) indicated that 13 of 22 genes 
were differentially expressed in at least one comparison 
between summed libraries. Results from the statistical 
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Figure 3 Overview of tiie salmon notocliord transcriptome from RNA-seq. (A) Genomic location of mapped reads in percentage in all 
salmon notochord libraries. (B) Distribution of Augustus-predicted gene length in basepair (bp) < 5000 bp; the others were not plotted. The red 
dotted line indicates the average gene length of all genes, which was 872 bp. (C) Distribution of expression (in log2 scale) of all genes in libraries, 
grouped by development stages (Tl:510±20 d°, T2: 610 ±20 d°, T3: 710 ±20 d°). (D) Venn diagrams showing number of differentially expressed 
genes from NOISeq between pair-wise comparison between summed spatial libraries (SumA, SumM, SumP) and summed temporal libraries 
(SumTl, SumT2, SumT3). Abbreviations: A (anterior), M (mid), P (posterior). 



analysis of differential gene expression in both methods 
are shown in Additional file 3: Table S2. Twelve of 13 
DEGs were confirmed by qPCR assay (except fos, which 
showed large variation between libraries). Among the DEGs 
confirmed with both methods, four DEGs {bmp6, hoxbSab, 
hoxclOaa, hoxdSaa) were significantly up-regulated from 
anterior to posterior (sumA-P) based on RNA-seq data. 
Four DEGs {tgfbly bmp6, dkk2, runxS) were significantly 
up-regulated from Tl to T3, while four DEGs ifetub, bmp4, 
col2al, coll la 1) were significantly down- regulated from 
Tl to T3. 

Expression of Hox genes within the notochord 

Pasta sequences of hox genes were further blasted 
against the NCBI nucleotide database, and annotated 
with full-length salmon hox genes with accession 
number provided in [15]. The expression of genes with 



the same annotation was fused before log transformation. 
In this study, 71 out of 118 hox genes in 13 clusters were 
expressed in the notochord (Figure 5A). Most of Hox 
clusters displayed the expected collinear expression, with 
some, such as members of the Hox A and B clusters, 
having the highest expression in the anterior samples, while 
others had the opposite pattern of expression, many of 
which were found in the C and D clusters (Figure 5A-C). 
The cluster with highest gene expression activity was the B 
cluster, while less active genes were found in A, C and D 
clusters (Figure 5A). Clustering of all hox genes revealed 
that about half of the 71 display a very low expression 
profile (Gl and G2 nodes in Figure 5C). Interestingly, 
several members of classical anterior clusters (A and B) 
contain many hox genes which have a posterior expression 
pattern (G5 node, in Figure 4C). Only a few genes had the 
highest expression in the anterior part of the notochord. 
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(See figure on previous page.) 

Figure 4 Validation of RNA-seq by qPCR. RNA-seq expression of 22 selected genes was confirmed by independent quantitative RT-PCR (qPCR) 
metliod. Expression pattern of four genes {fetub, tgfbl, bmp4, bmp6) involved in transforming growth factor-beta (TGF|3)/bone morphogenic protein 
(BMP) signaling pathways (A); four genes {dkk2, fos, t7llo, Irp6) in Wnt signaling pathway (B); three hox genes {hoxbdob, hoxlOoo, hoxdSoo) (C); 3 
Hedgehog genes {shh, ihhb, twhh) (D); three transcriptional factors {runx3, atf4, sox9) involved in osteochondro-progenitors lineage (E); three fibrillar 
collagen genes (col2a1, colllal, col11a2) (F); two calcium binding proteins (sprc, si 00a) (G) across all nine notochord libraries (T1 A, Tl M, Tl P, T2A, 
T2M, T2P, T3A, T3M, T3P) at three developmental stages (Tl :510 ± 20 d°, T2: 610 ± 20 d° T3: 710 ± 20 d°). Of the expression of individual genes in each 
panel, the X-axis indicates notochord libraries; the Y-axis indicates the gene expression in one library as fold change relative to Tl A library (Tl A: anterior 
part of salmon notochord at 510 + 20 d° stage), measured with RNA-seq (cyan) and qPCR (red). Error bar indicates mean ± SD. Overall, Abbreviations: A 
(anterior), M (mid), P (posterior). Gene name: Fetuin-B (fetub); Transforming growth factor beta-1 {tgfbl); Bone morphogenetic protein 4 {bmp4); Bone 
morphogenetic protein 6 {bmp6); Low-density lipoprotein receptor- related protein 6 (/rp6); Dickkopf-related protein 2 (d/c/c2); Transcription factor 7-like 
1-A {tyilo); Proto-oncogene c-Fos (fos); Homeobox protein HoxB6ab (hoxbdob); Homeobox protein HoxClOaa (hoxclOoo); Homeobox protein HoxD3aa 
(hoxdSoo); Indian hedgehog B protein (ihhb); Sonic hedgehog protein A (s/7/7);Tiggy-winkle hedgehog protein (tw/ih); Transcription factor SOX-9 (sox9); 
Runt-related transcription factor 3 {runxS); Cyclic AMP-dependent transcription factor ATF-4 (c/tf4); Protein SlOO-Al {si 00a): SPARC(sprc); Collagen 
alpha-l(ll) chain (co/2a/); Collagen alpha-l(XI) chain (co// /a/); Collagen alpha-2(XI) chain {colllo2). 



including; hoxbSba, hoxc4bb, hoxb4bb, hoxa4aa, hoxb2aa, 
hoxbSab and hoxbSba (G4 node, in Figure 5C). Also, in 
the whole dataset, there seems to be an overrepresen- 
tation of hox genes having a posterior expression pattern 
(Figure 4C). For most of the genes, expression was 
consistent throughout all three developmental stages. 
However, in the highly expressed hox genes (G5 node, 
in Figure 4C), gene expression increased at 610 d°. 
The posterior expression of a 3' located hox gene 
(hoxdSaa), hoxb6ab and another 5' located hox gene 
(hoxclOaa), were validated by qPCR (Figure 4C). 

Expression of key molecules in chondrogenic lineage, 
ECM formation and mineralisation 

To determine if there exist common pathways that 
determine the differentiation of osteoblast, chondroblast 
and chordoblast (notochord epithelium cells), and thus a 
possible evolutionary link in tissue morphogenesis, key 
genes directing mesenchyme stem cells toward osteogenic 
and chondrogenic lineages were searched, and the expressed 
genes are listed in Figure 6. Runx2 {cbfal), osterix and 
osteocalcin {bgp), which are genes known to be involved in 
osteogenesis, were not detected in our analysis, and no 
expression of runx2 was confirmed by qPCR. Interestingly, 
runxl and runxS, two genes closely related to runx2, but 
involved in the chondrogenic lineage, were expressed 
{runxS expression is shown in Figure 6). The master 
regulator of chondrogenesis, sox9, was highly expressed in 
the notochord, confirmed by qPCR (Figure 4E). Transcripts 
of cofactors such as soxS, sox6 that bind to sox9 in order to 
activate col2al transcription, were found (Figure 6). In terms 
of chondrogenic cell differentiation, transcripts of many 
factors {sox9, sox6, soxS, ihhb, tgfbS, col2al) were detected 
(Figure 6). A strong up-regulation of tgfbl, confirmed with 
qPCR, was observed in T3 (Figure 4A). 

Overall activity in the notochord during ontogeny 

KEGG analysis of differentially expressed pathways between 
sumTl-T2 (Figure 7A) and sumT2-T3 (Figure 7B) 



revealed that most pathways were down-regulated 
between T2 and T3. Considering pathways with 
above 10 genes being regulated (marked by dotted line in 
Figure 7), comparisons revealed for up-regulated genes 
T1-T2: osteoclast differentiation, MAPK signalling 
pathway and purine metabolism; T2-T3: Neuroactive 
ligand receptor interaction. The numbers of up-regulated 
genes were few compared to those down-regulated, 
while the analysis revealed no pathways for T1-T2, 
and a comparison of T2 and T3 showed many pathways 
being down-regulated, including focal adhesion, ECM 
receptor interaction, cell cycle, protein digestion and 
absorption, protein processing in the endoplasmic 
reticulum etc. (Figure 7B). The canonical Wnt pathway 
was down-regulated from the T2 to T3 stage (Figure 7B), 
and the full map with DEGs plotted is shown in 
Additional file 4: Figure S3A. Interestingly, two genes 
(KEGG ID) repressing the canonical Wnt pathway 
{sfrpl and dkk2) and another gene {smad4), were 
up-regulated from T2-T3, and 11 other genes {wnt2, 
lrpS/6, smad4y EP300, fosll) were down-regulated, 
suggesting a complete picture of down-regulation of 
the pathway from T2-T3. Another KEGG pathway 
(protein processing in the endoplasmic reticulum) consisted 
of 31 down-regulated DEGs and only two up-regulated 
genes, as shown in Additional file 4: Figure S2B. 

Spatial expression of co/7 1a2 in the notochord 

The notochord in the stages studied is composed of a 
collagenous sheath, a single layer of chordoblasts, and a 
central core of vacuolated cells - the chordocytes. At 
500 d° (Figure 8A), the notochord is morphologically 
unsegmented, while at 700 d°, early segmentation becomes 
visible (Figure 8B). The changes in the transcriptome that 
occurs between 500 and 700 d° may largely represent a 
preparation for the profound structural changes that take 
place, particularly in the sheath when mineralisation com- 
mences, ultimately to form the chordacentrae (Figure 8E). 
The sheath increases substantially in the intervertebral 
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Figure 5 Hox gene expression from RNA-seq within the salmon notochord during segmentation. 71 out of 1 18 salmon iiox genes were 
expressed and tine average expression (on log2 scale) of hox genes in all salmon notochord libraries (A). (B) Heatmap showing RNA-seq expression 
(in log2 scale) of hox genes in nine notochord libraries; x-axis primarily grouped by developmental stage (Jl:510 ± 20 d° T2: 610 ± 20 d° T3: 710 ± 20 d°); 
y-axis grouped by four major clades in hox clusters. (C) K-mean clustered heatmap showing RNA-seq expression (in log2 scale) of hox genes in 
nine notochord libraries. Lowly expressed posterior genes were represented in Gl; generally lowly expressed genes were grouped in G2; most 
highly-expressed genes were grouped in G3; G4 represented anteriorly-expressed genes and G5 shown the highly-expressed hox genes either in 
mid- and posterior segments or had no clear pattern. The horizontal clustering indicates a closer distance in overall expression pattern of hox genes 
between temporal libraries (Tl, T2, T3) and, to a further distance between spatial libraries (A, M, P). Abbreviations: A (anterior), M (mid), P (posterior). Note 
that a total of 1 10 hox genes were plotted in all figures (A-C), since the rest 8 hox genes are pseudogenes. 
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Figure 6 Clustered heatmap of chondrogenic factors expressed in the notochord. Gene networks that possibly direct chondrogenic 
differentiation pathway from mesenchyme stem cells were expressed in the notochord RNA-seq dataset. Genes were filtered with a threshold of 
blast score > 300 and represented by genes with highest blast score if they had the same gene symbol. Gene expression from RNA-seq from three 
spatial segments were summed and represented as sumTl, sumT2 and sumT3 (Tl:510±20 d°, T2: 610 ±20 d°, T3: 710 ±20 d°) and plotted as fold 
change to the average among the summed libraries. Genes were clustered according to their expression pattern. 
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Figure 7 Actively regulated KEGG pathways in salmon notochord during segmentation. Numbers of differentially expressed genes (DEGs) 
were grouped by KEGG ID in pathways. Only genes with a blast score over 300 and fold change >1.5 from NOISeq were counted as DEGs and plotted in 
KEGG. Up-/down-regulated KEGG pathways from sumTl-T2 (510 ±20 d°-610±20 d°, pre-seg mentation) were listed in (A). Up-/d own-regulated KEGG 
pathways from sumT2-T3 (610 ± 20 d°- 710 ± 20 d°, segmentation) were listed in (B). Up-regulated pathways were ranked by the ratio of the number of 
up-/down-regulated genes (shown if the ratio is above 1), and vice versa for down-regulated pathways. Pathways with fewer than 10 DEGs (dotted 
vertical line) were considered to be insufficient for inclusion. Full maps of two pathways (asterisk marked) are illustrated in Additional file 4: Figure S3. 



regions (Figure 8E-F), forming a thick intervertebral 
ligament, while in the region of the chordacentrae, 
growth in thickness of the sheath ceases (Figure 8F), 
In situ hybridization was performed on longitudinal 
cross-sections of notochord at 710 d° and spatial 



expression of collla2 was found segmentally in the 
chordoblast layer, while no expression was found in 
the chordocytes (Figure 8C). No background signal was 
found using the sense probe of the gene (Figure 8D), nor 
was segmental expression of collla2 observed at 510 d° 
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Figure 8 Notochord morphology and in situ hybridization of col11a2. (A, B) show the morphology before, during and after mineralisation 
start in the notochord sheath. (A) The non-segmented notochord at 500 d° Arrow: notochord sheath; arrowhead: chordoblasts; ca: cartilaginous arcualia. 
Scale bar: 100 um. (B) Notochord at start of mineralisation in the sheath (arrow), and slightly segmented chordoblasts (arrowhead). Scale bar: 100 um. 
(C) Tissue-specific, segmental gene expression of coll la2 antisense in chordoblasts (arrowhead) in the notochord at 710 d°. No staining in the 
notochordal sheath (arrow). Scale bar: 20 |jm. (D) No background signal was found using the sense probe of coll lo2 (710 d°). Arrowhead: chordoblasts, 
arrow: notochord sheath. Scale bar: 20 pm. (E) Notochord and surroundings, 1000 d°. Red arrow indicates the mineralised chordacentrum, the earliest 
vertebral element. Ca: arcualial cartilage. (F) Differentiated chordacentrum (red arrow), 1300 d°. Notice the thickened notochordal sheath between 
segments (il: intervertebral ligament), and the thin remaining part of the sheath. Ca: arcualial cartilage. Scale bar: 100 |jm. 



(data not shown). This suggests differentiation of a distinct 
population of cells that expressed collla2 during 
notochord segmentation. Both RNA-seq and qPCR 
have confirmed that colllal was down-regulated from 
T2 to T3, although inconsistent results were obtained 
between the two methods for collla2. 

Building and mineralisation of the notochord sheath 

To further elucidate which proteins build up the 
notochordal sheath, extracellular matrix (ECM) genes 
were retrieved from the annotated RNA-seq pool using 
Gene Ontology. Only genes with more than 100 reads in 
the highest coverage of all libraries were further clustered 
(k-means) over time (T1-T3), in order to find patterns in 
expression of ECM proteins for the pre-segmented 
notochord (510-610 day"), and during segmentation 



of the notochord (610-710 day"). K-mean clusters 
were filtered, and DEGs confirmed by NOISeq were 
demonstrated in all ECM genes. In genes that were 
up-regulated from Tl to T3 only six annotations were 
identified {ctgfl, acm2, efemp2, fbln7, tfp2 and tgfbl, 
Figure 9A). In genes down-regulated from Tl to T3, 
as many as ten with the same annotation were found; 
bean, colllal, col24al, col2al, crtap,Jmod,fnl, leprel, lox 
and loxl (Figure 9B). 

Discussion 

About 65000 salmon gene transcripts were expressed in 
the notochord during the developmental stages analysed. 
This is an expected number of transcripts, considering 
the partial tetraploidy of salmonids [29]. In for example 
zebrafish, which have undergone only one teleost-specific 
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Figure 9 Differentially expressed ECM genes in the salmon notochord. Genes related to ECM annotated from Gene Ontology (GO) were merged 


with differentially expressed gene list output from NOISeq with threshold 0.8 in summed temporal libraries T1-T3 comparison, thereby showing differentially 


expressed ECM genes from Tl to T3 in (A,B). Heatmap of RNA-seq expression of these genes in SumTl, SumT2 and SumT3 notochord libraries were 


showed as fold change to the average expression of all summed libraries. Up-regulated ECM genes were listed in (A) and down-regulated genes in (B). 



genome duplication in the course of evolution, -51000 
transcripts were identified in a tissue gene expressional 
study using a similar genome-mapping method [30]. De 
novo assembly of the RNA-seq data could have been 
implemented, but this method may produce a much 
higher number of predicted transcripts, as has been 
observed in rainbow trout {Oncorhynchus my kiss) and 
catfish {Ictalurus punctatus), where de novo assembly 
produced about --130000 and -370000 predicted 
transcripts respectively [31,32]. With an average predicted 
coding sequence length of around 872 bp, the results are 
comparable to what was predicted using de novo assembly 
methods in catfish [32], but much longer than that 
obtained in rainbow trout [31]. Of the -65000 predicted 
transcripts obtained in the present study, annotation was 
retrieved for -55000 transcripts (83%), which is a much 
higher rate compared to rainbow trout, where only 50% 
transcripts were annotated [32]. However, due to the 
fragmented structure of the current salmon genome 
assembly, which consists of a large Sanger-sequenced 



data-set with a contig N50 of 9342 bp, it was difficult 
to predict whole protein lengths and therefore also 
splice variants in a large fraction of the expressed genes. 
Moreover, due to the recent genome duplication in 
salmonids, many paralogous genes display high sequence 
similarity, which may cause cross-mapping, leading to 
assembly errors. 

Quantitative PGR was used to validate gene expression 
profiles. Twelve out of thirteen DEGs from NOISeq were 
confirmed by qPGR, which indicates the high reliability of 
NOISeq for detecting DEGs. qPGR appeared to be more 
sensitive than RNA-seq in detecting differentially expressed 
genes, which was expected due to presence of isoforms 
which share mapping unequally, a fragmented genome and 
lack of replicate samples in RNA-seq [33]. Interestingly, 
there was a high correlation in gene expression data among 
the hox clusters, as they all displayed collinearity in their 
expression patterns, in spite of some of them being 
expressed at a very low level. This gene family consist of 
very well-characterized full-length genes identified by 
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manual cloning [15]; therefore a fragmented genome 
might not have affected gene expression profiles of 
these genes. The lack of replicates also made it difficult to 
calculate differential gene expression, but by using NOISeq 
we were able to calculate differential gene expression 
profiles between stages by combining the expression values 
of three spatial libraries (A, M, P). NOISeq was chosen 
mainly because it can simulate technical replicates when 
few or no replicates are available, which is the case of this 
study [21]. Since NOISeq does not rely on parametric 
assumptions, it is also more effective in both identifying 
differentially expressed transcripts with low read counts 
and reducing the rate of false discoveries [21]. 

Of all the 118 salmon hox genes including 8 pseudogenes, 
71 were expressed in the developing notochord, and also 
coUinear expression was observed in corresponding hox 
clusters along the anterior-posterior axis. Interestingly, our 
study shows that duplicated Hox cluster pairs in salmon 
have different gene expression patterns, and also this is 
most obvious in HoxBa clusters (Figure 5A). This fact could 
be linked slow evolution of HoxBa, both in number of 
genes and in their gene and intergenic sequences [15]. Our 
data may therefore suggest that HoxBa have evolved 
different functional roles in the notochord. So far, only one 
study in teleosts have previously identified expression of 
four hox genes in the notochord, among which hoxbS dis- 
played the most posterior and highest expression within the 
notochord [14]. In the notochord, we also found the highest 
expression in HoxB clusters. Another RNA-seq study of 
whole eel embryos {Anguilla Anguilla), which includes the 
notochord, has also shown that the HoxBa has the highest 
expression [34]. Furthermore, clustering of all hox expres- 
sion patterns in the notochord indicated that highly- 
expressed hox genes increase in expression at 610 day, and 
this expression is retained until 710 day (G5 node, in 
Figure 5C). Interestingly, the expression levels of all hox 
genes were more similar in the medial part of the noto- 
chord throughout the study period than in other regions, 
which is in accordance with the observation that the initial 
chordacentrae develop there [9,12,13]. Another point was 
that several genes found in the A and B clusters show pos- 
terior expression patterns. Prince et al. [14] have observed a 
similar posteriorisation of hox gene expression in the noto- 
chord compared to surrounding tissues, including the CNS 
and paraxial mesoderm. Here, there are two plausible expla- 
nations; either hox genes have a more posterior fate in the 
notochord, or the posterior end of the notochord represents 
its own patterning unit. In the latter case, we would not be 
able to detect collinearity in the posterior notochord, since 
all this tissue was in the posterior sample. But as most of 
the hox genes display a posteriorisation of expression, it is 
likely that the latter assumption is to some extent correct. 

A wide search for factors implicated in tissue mineralisa- 
tion resulted in finding of chondroblast-specific "master 



genes" expressed in notochord such as sox6, soxS and sox9 
[35], in addition to other genes associated with the fate of 
chondroblasts, such as tgfbS, ihhb and col2al. However, no 
genes associated with the fate of specific osteoblasts, such 
as runx2, osterix and bgp^ were found. The chondrogenic 
genes SoxS and Sox6 have previously been shown to be 
expressed in the zebrafish notochord. Knock-out of both 
genes in zebrafish inhibits ECM formation in the 
notochord, which subsequently induces apoptosis in 
notochord cells [36]. Interestingly, the sox9a/b mutant 
in zebrafish displays a truncated notochord, supporting 
the notion of a function of sox9 in chordocytes [37]. Ihhb 
is essential for chondrogenesis in vertebrates, including 
teleost fish, where its expression has been associated with 
chondrogenic tissue [38]. In the same study, ihhb was 
identified in the notochord, further supporting the concept 
of evolutionary closeness between chondroblasts and chor- 
doblasts. TgfbS was another gene found to be expressed in 
the notochord, and this is also expressed in a number of 
cartilage tissues in zebrafish [39]. TgfbS is also highly 
expressed in the notochord during early development in 
the zebrafish, but as development progresses it ceases to be 
expressed in the notochord. Similarly, we observed a 
reduced expression of tgfbS over time in salmon notochord. 
One of the most highly expressed genes within the 
notochord is the major collagen of hyaline cartilage, col2al 
[40]. This collagen type displays a much higher level of 
expression than the major type I collagen {coll) found in 
bone, further confirming the chondroblast-like nature of 
chordoblasts [13]. Since the notochord appeared early in 
chordate evolution, as a key structure characterising the 
clade, and anticipating by far the emergence of hyaline 
cartilage, which probably first appeared in jawed fish, it 
may be fair to say that chondroblasts are chordoblast-like, 
having co-opted morphogenic pathways that evolved 
together with the notochord in the pre-chordate lineage, 
early in the Cambrian (see discussion in [12]). 

Both qPCR and RNA-seq employed in this study show a 
significant down-regulation of the fibril collagen colllal at 
onset of mineralisation (610d° -710d°), suggesting an overall 
reduced synthesis and re-modelling of extracellular matrix 
in the notochord sheath concurrent with segmentation and 
mineral crystal deposition. Nine other ECM genes were also 
down-regulated at onset of mineralisation (610d° -710d°); 
these included bean, col24a\, col2al, crtap, fmod, fnl, 
leprel, lox and loxL It was also shown that at 710d°, 
collla2 was spatially expressed in regions not undergoing 
mineralisation - the sites of the prospective intervertebral 
joints - which suggests that the overall down-regulation 
of genes may be due to compartmentalisation of 
expression (Figure 8C). Likewise, expression of ALP 
in chordoblasts in segments where mineralisation occurs 
substantiates such a notion. In the notochord, it has been 
observed that induction of ECM-specific genes over time 
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(ctgf, efemp2 or fibulin 4, fbln? and tfpi2), and proteins 
encoded by these genes, might further contribute to the 
specificity of the segmentation and mineraUsation 
processes. In the KEGG analysis of differentially expressed 
genes performed in the present study, we discovered a 
massive down-regulation of many pathways between 
510-710 day". However, only the Wnt pathway was 
found to be down-regulated early, in the period between 
610 and 710 day °, which may indicate that the Wnt 
pathway is involved in the segmentation and subsequent 
mineralisation of the chordacentrae (Figure 7B). Indeed, 
key Wnt signaling pathway genes were affected, including 
fzd2, fosU nrlc2, wnt2, lrp5/6, fosU sfrpl and dkk2 
(Additional file 4: Figure S2A). 

Experimental manipulations in the chick embryo have 
excluded the possibility that the notochord plays a role 
in segmentation of the vertebral column, indicating that 
subdivision of the sclerotome alone leads to this patterning 
in avian species, and probably in all amniotes [41]. A 
study performed on the Japanese medaka {Oryzias latipes), 
supports the notion that this mechanism is also involved 
in teleost species [42]. Prior to mineralisation of the 
notochord sheath in medaka, sclerotome-derived cells, 
which express collOal, migrate to surround the notochord 
at segments where the prospective chordacentrae are to 
develop. The cells differentiate into osteoblasts, and it has 
been suggested that they, rather than the notochord 
epithelium, secrete molecules that induce mineralisation 
of the sheath, thus forming the chordacentrae and 
subsequently depositing osteoid to form perichordal 
vertebral bone. Such a mechanism, in which segmentation 
is independent of the notochord, would be similar to 
that described in amniotes [41]. On the other hand, 
in zebrafish, sclerotomal cells do not migrate to surround 
the notochord in strictly segmented compartments [43], 
and, in salmon, the notochord epithelium expresses 
ALP - a key enzyme in tissue mineralisation - in a 
spatio-temporal manner that coincides with development 
of the chordacentrae; here osteoblast precursor cells 
surrounding the notochord do not express ALP, and thus 
do not appear to play a role in the mineralisation of the 
notochord sheath. In medaka, an alternative mechanism 
for the accumulation of osteoblast precursors at specific 
sites along the notochord [42] may be that the cells 
respond to unidentified chemoattractant guiding factors 
that are segmentaly secreted by the notochord epithelium. 
Indeed, to ensure an uninterrupted development of 
the vertebrae, these cells need to be in place before the 
chordacentrae form. Moreover, in zebrafish, when the 
osteoblasts surrounding the notochord are abolished, the 
initial development of chordacentrae proceeds normally, 
but they continue to grow, so that the spaces separating 
adjacent centrae become reduced, in some cases finally 
leading to partial fusion [8]. Hence, osteoblasts that are to 



form perichordal bone may provide signalling factors 
that inhibit the notochord epithelium in its secretion 
of molecules that promote mineralisation of the sheath, 
arresting further growth of the chordacentrae, and, by 
extension, defining a boundry with non-mineralised 
segments that correspond to where prospective interverte- 
bral ligaments are to develop. Here, we have identified a 
metameric expression of collla2 within the notochord that 
colocalizes with these segments. Apart from being an 
extracellular structural matrix protein found in the 
core of the collagen type II fibril, collagen type XI 
may inhibit maturation of osteoblasts [44] in the region 
where intervertebral ligaments develop - here, instead of 
osteoblasts, fibroblasts condense around the notochord to 
form collagenous tissue [11]. One of the end players in the 
Wnt pathway, which is down-regulated at the onset of 
segmentation, is collagen type XI. This molecule may act 
as a downstream inhibitor of mineralisation along the 
whole notochord prior to segmentation, as expression of 
collla2 is, at the onset of segmentation, confined to 
regions where mineralisation does not occur. Hence, 
through controlling inhibitory signals, the Wnt pathway 
may play a key role in notochord segmentation. 

The gathering of osteoblast precursors around the 
notochord, and the inhibitory effect of osteoblasts on 
growth of the chordacentra, suggests a reciprocal 
interaction between notochord and surrounding tissue. 
However, our finding of a complex dynamic regulation of 
a range of anterior-posterior polarity genes and the 
segmental expression of collla2 within the notochord 
epithelium supports the notion that autonomous segmen- 
tal patterning processes exist in the teleost notochord. 

Conclusions 

Sequencing and partial analysis of the notochord 
transcriptome of salmon, the first dataset of its kind 
from a vertebrate, has provided new insight into the 
cellular and molecular genetic mechanisms of this 
organ. The coUinear expressions of 18 hox clusters in 
the notochord represented by 71 hox genes, suggest a 
possible role of hox genes in anterior-posterior patterning 
of the notochord. When notochord segmentation occurs, 
we find an extensive down-regulation of genes, which may 
reflect a shutdown of complex regulatory developmental 
pathways leading up to this event. Expression of collla2 
was detected in segments where prospective intervertebral 
ligaments develop, and we hypothesise that the Wnt 
pathway, where collagen type XI is one of the end players, 
may play a key role in segmental patterning. Furthermore, 
the expression of genes that determine differentiation 
of chondroblasts indicates that cartilage may have 
evolved after the co-option of morphogenic pathways 
that initially evolved together with the notochord in the 
pre-chordate lineage. 
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Additional files 



Additional file 1: Table SI. RNA library and reads mapping summary. 
Quality of RNA libraries and reads, as well as mapping statistics 
summarised in one table. 

Additional file 2: Figure SI. Quality check summary. (A) Quality score 
of nucleotide positions in short reads from a titration run. (B) Nucleotide 
distribution in short reads from a titration run. 

Additional file 3: Table S2. Primer sequence and additional 
information of selected genes in qPCR validation and in situ 
hybridization. Detailed information about 23 selected genes used in qPCR 
validation of RNA-seq gene expression. One-way ANOVA was used to 
define significance of differentially expressed gene (DEG) if p<0.05 among 
summed libraries (both spatially and temporally). DEG analysis using 
NOISeq was used to detect genes from RNA-seq dataset. 

Additional file 4: Figure S2. Full KEGG pathway maps for two 
important pathways. Wnt signaling pathway and protein processing in 
endoplasmic reticulum are illustrated in (A) and (B), respectively. 
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